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1 Introduction 

In this report, we present numerical experiments with three particular ABS algorithms: 

(1) The Huang or modified Huang algorithm, 

(2) The implicit LU algorithm, 

These algorithms have been used for finding a solution of the linear KKT system 
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where A G i?™'", B G i?"'" is symmetric, b G K"-, c G i?™, x G K^, y G and m < n. 

The considered algorithms belong to the basic ABS class of algorithms for solving a linear system 

Ax = b, which is a special case of the scaled ABS class, given by the following scheme: 

(A) Let xi G i?" be arbitrary and Hi G i?"'" be nonsingular arbitrary. Set i = 1. 

(B) Compute the residual = Axi — b. If r, = 0, then stop, solves the problem. Otherwise 
compute Si = HiA^Vi, where Vi G i?" is arbitrary, save that vi, . . . ,Vi are linearly inde- 
pendent. If Si 7^ 0, then go to (C). If Sj = and rfvi = 0, then set Xj+i = Xi, -fTj+i = Hi 
and go to (F). If Sj = and rfvi ^ 0, then stop, the system is incompatible. 



(C) Compute the search vector pi by 



Pi 



where Zi G i?" is arbitrary, save that zjHiA^Vi / 0. 
(D) Update the estimate of the solution by 

Xij^l ^i ^iPi-i 

where the stepsize Oj is given by 

cui = rfvi/pfA^Vi. 
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(E) Update the Abaffian matrix by 



Hi+i = Hi- HiA^v^wfHi/wfHiA^Vi, 
where Wi G i?" is arbitrary, save that wfHiA^Vi ^ 0. 

(F) If i = m, then stop, Xj+i solves the problem. Otherwise increment the index i by one and 
go to (B). 

From (E), it follows by induction that Hi^iA^Vi = and Hf^^Wi = 0, where Vi = [vi, . . . , Vi] 
and Wi = [wi, . . . ,Wi]. One can show that the implicit factorization V^APi = Li holds, where 
Pi = bi;--- }Pi] s-iid Li is nonsingular lower triangular. Moreover the general solution of the 
scaled subsystem Ax = V^b can be expressed in the form 

x = Xi + Hfq, (1) 

where q E i?" is arbitrary (see ||2| for the proof). 

The basic ABS class is the subclass of the scaled ABS class obtained by taking Vi = e^, 
being the i-th unitary vector (i-th column of the unit matrix). In this case, residual rj = Axi — b 
need not be computed in (B), i-th element rfci = ajxi — hi suffices. 

This report is organized as follows. In Section 2, a short description of individual algorithms 
is given. Section 3 contains some details concerning test matrices and numerical experiments. 
Section 4 discusses the results of the numerical experiments. The Appendix contains all numer- 
ical results. For a listing of the ABS source codes see [12]. For other numerical results on ABS 
methods see [1,11]. 



2 ABS algorithms for KKT equations 

In this report, we use three linear ABS algorithms, two belonging to the basic ABS class and one 
belonging to the scaled ABS class, for solving the KKT linear systems exploiting their structure. 
To simplify description, we will assume that A € K^''^, m < n, has full row rank so that Si ^ 
in (B). First we describe the used linear solvers. 

The Huang algorithm is obtained by the parameter choices Hi = I , Vi = Ci, Zi = ai, Wi = ai. 
Therefore 

Pi = HiUi (2) 

and 

H,+i = Hi- Pipf/ajpi, (3) 
From (2) and (3), it follows by induction that pi E Range{Ai) and 

H,+i = I-P,Dr^P^. (4) 

where Ai = [ai,...,aj]. Pi = [pi,...^pi] and Di = diag{aiPi, . . . ,afpi). Moreover Hi is 
symmetric, positive semidefinite and idempotent (it is the orthogonal projection matrix into 
Null{Ai-i)). Since the requirement pi G Null{Ai-i) is crucial, we can improve orthogonality by 
iterative refinement pi = Hjai, j > 1 (usually j = 2), obtaining the modified Huang algorithm. 

The Huang algorithm can be used for finding the minimum-norm solution to the compatible 
underdetermined system Ax = 6, i.e. for minimizing s.t. Ax = b. To see this, we use the 
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Lagrangian function and convexity of Then a; is a required solution if and only if a: = A^u 
for some u G or x G Range{A^). But 

m 

Xm+1 =Xl-^ aiPi 
i=l 

by (D) and pi G Range{Ai) C Range{A^), so that if xi = 0, then Xm+i £ Range{A^). 

A short description of two versions of the Huang and modified Huang algorithms follows. 

Algorithm 1 

(Huang and modified Huang, formula (3). 
Set xi = and Hi = I. 
For z = 1 to m do 

Set Pi = Hitti (Huang) or 

Pi = Hi^HiQi) (modified Huang), 

di = afpi and x^+i = Xj - {{afxi - hi)/di)pi. 

If ? < m, then set -ffj+i = Hi— Pipf /di. 

end do 



Algorithm 2 

(Huang and modified Huang, formula (4). 
Set xi = and Pq empty. 
For i = 1 to m do 

Set Pi = {I- Pi-iD-\P^_^)ai (Huang) or 

Pi = {I- Pi-iD-\Pl_^){I - Pi_iD-\Pl^)ai (modified Huang), 
di = ajpi and Xj+i = Xi - {{afxi - bi)/di)pi. 
If i < m, then set Pj = {Pi-i,pi]. 
end do 



The implicit LU algorithm is obtained by the parameter choices Hi = I, Vi = Ci, Zi = Ci, 
Wi = Ci- Since W[f Hi j^i = [Ij, 0]"^i?i_|_i = 0, the first i rows of the Abaffian matrix -ff^+i must 
be zero. More precisely, the Abaffian matrix has the following structure, with Ki G RJ^' 



-1,1 



H 





Ki Ir. 







(5) 



Only the matrix Ki has to be updated. Expression (5) implies that only the first i components 
of the vector pi = Hjci can be nonzero and the i-th component is unity. Hence the matrix Pi is 
unit upper triangular so that the implicit factorization A = LP~^ is of the LU type with units 
on the diagonal. 

A short description of two versions of the implicit LU algorithm follows (the second version 
is also called the implicit LX algorithm). 

Algorithm 3 

(Implicit LU with explicit column interchanges). 
Set xi = and Hi = I. 
For z = 1 to m do 

Set Si = Hitti 

(only (i — l)(n — i + 1) nonzero elements of Hi are used). 
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Determine di = = m.aj<.{\sfej\, j = 1, . . . ,n). 
(only n — i + 1 nonzero elements of Sj arc used). 

If ki ^ i, then swap columns of A and elements of x and s with these indices. 

Set Xi+i = Xi- {{ajxi - bi)/di)Hfei 

(only i nonzero elements of Xj are updated). 

If i < m, then set -ffj+i = Hi — SicJ HJ / di 

(only i{n — i) nonzero elements of -ffj+i are updated). 

end do 

Algorithm 4 

(Implicit LX algorithm, or implicit LU with implicit column interchanges). 
Set xi = 0, Hi = I. 
For z = 1 to m do 

Set Si = HiQi 

(only (i — l)(n — i + 1) nonzero elements of Hi are used). 

Determine dfc. = = max.{\sfej\,j = 1, . . . ,n). 

(only n — i + 1 nonzero elements of s j are used) . 

Set Xi+i = Xi- {{aJxi - h) / dki)Hj 

(only i nonzero elements of updated). 

If i < m, then set -ffj+i = Hi — Sie^,Hj / dk^ 

(only i{n — i) nonzero elements of -ffj+i are updated). 

end do 

We now consider how to use the above algorithm for solving KKT systems, exploiting their 
structure. Consider the linear KKT system 

Bx + A^y = b, (6) 

Ax = c. (7) 

The underdetermincd equation (7) can be solved by an arbitrary ABS method obtaining the 
particular solution x^+i- Since H^+iA^ = 0, we get 

Hm+iBx = Hm+ib, (8) 

by multiplying (6) by Hm+i- The coupled system (7)- (8) is overdetermined but compatible so 
that its unique solution can be found by an arbitrary ABS method. Since rank{Hjn+i) = n — m, 
m equations will be eliminated. Notice that solution of (7)- (8) is obtained in two steps. First 
we solve (7) to obtain Xm+i and Hm+i- Then we construct Hm+iB and H^+ib and starting 
with Xm+i and Hm+i we continue with the ABS method to solve (8). 

If we use the implicit LU algorithm, then matrix H^+i has the special form (3) and the first 
m equations of (8) are trivially satisfied, leaving a determined system consisting of (7) and the 
equation 

SrnBx = Smb, (9) 

where Sm = [Km, In-m]- Moreover, substituting the general solution, see (1) 



with q G arbitrary into (6) and using the special form of Hm+i, we obtain 





SjnBS^ 











(12 







Sm{b - BXm+l) 



(10) 
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Thus qi can be chosen arbitrarily, e.g. qi = 0, and q2 is a solution to the n — m dimensional 
system 

SmBS'^q2 = Sm{b - Bxm+i), (11) 

which can be solved by any ABS method. 

Once X is determined, we obtain y by solving the compatible overdetermined system A^y = 
b — Bx. Since the equation Ax = c was solved beforehand, we are in the same position as in the 
case of the least-squares solution of an overdetermined linear system. Therefore, using P^, we 
can construct the lower triangular matrix = AP^n and get the solution from the equation 
Lly = Pl{b-Bx). 

3 Description of computational experiments 

Performance of ABS algorithms has been tested by using several types of ill-conditioned matrices. 
These matrices can be classified in the following way. The first letter distinguishes matrices with 
integer T and real 'R' elements, both actually stored as reals in double precision arithmetic. The 
second letter denotes randomly generated matrices 'R' or matrices determined by an explicit 
formula 'D'. For randomly generated matrices, a number specifying the interval for the random 
number generator follows, while the name of matrices determined by the explicit formula contains 
the formula number (Fl, F2, F3). The last letter of the name denotes a way for obtaining 
ill-conditioned matrices: 'R' - matrices with nearly dependent rows, 'C - matrices with nearly 
dependent columns, 'S' - nearly singular symmetric matrices, 'B' - both matrices in KKT system 
ill-conditioned. More specifically: 

IR500 Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500]. 

IR500R Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500] perturbed in addition to have two rows nearly dependent. 

IR500C Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-500,500] perturbed in addition to have two columns nearly dependent. 

RRIOO Randomly generated matrices with real elements uniformly distributed in the interval. 
[-100,100] 

IDFl Matrices with elements aij = \i ~ j\, 1 < i < m, 1 < j < n (Micchelli-Fiedler matrix). 
IDF2 Matrices with elements Uij = \i — l<i<m, l<j<n. 
IDF3 Matrices with elements aij = \i + j — [m + n)/2|, 1, < i < m, 1 < j < n. 
IR50 Randomly generated matrices with integer elements uniformly distributed in the in- 
terval [-50,50]. 

Matrices with linearly dependent rows were obtained in the following way. The input data 
contain four integers which specify two row indices i\, i2, one column index and one exponent 
i4. Then the row is copied into ai,^. Furthermore Oj^jg is set to zero and 0^223 to 2"*". Similar 
procedures are used for columns and symmetric matrices. 

Right hand sides were determined from the given solution vectors x* by the formula b = 
Ax*. Solution vectors were usually generated randomly with integer or real elements uniformly 
distributed in the interval [-10,10]. 

We have tested KKT systems with m » n/2, m = n/2 and m « n/2. These systems 
were solved by using the modified Huang algorithm applied to the coupled system (7)- (8) and by 
two variants of the implicit LU algorithm with explicit column interchanges. The first variant is 
intended for solving the coupled system (7), (9) while the second one uses (7) together with the 
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transformed system (11). For comparison, we have tested three additional methods implemented 
by using the LAPACK routines. 

The first method, denoted by lu lapack, is in fact a direct solution of a complete KKT system 
by using the Bunch-Parlett decomposition. This is carried-out by the LAPACK routine DSPSV. 
Notice that we deal with an n + m dimensional indefinite system in this case. 

Another approach, a range-space method, is based on the Bunch-Parlctt decomposition of 
the (possibly indefinite) matrix B. Using the LAPACK routine DSPSV, we obtain the sohition 
of the matrix equation B[A^, b] = [A'^, b] together with the decomposition B = LDL^, where L 
is n-dimensional lower unit triangular and D is n-dimensional block diagonal with the blocksize 
1 or 2. Then the symmetric matrix C = AA^ = AB^^A^ is buih, and the sohition y to the 
equation Cy = 6 — c is found, again by using the LAPACK routine DSPSV. Finaly, we solve the 
system LDL^x = b- A^y by using the LAPACK routine DSPTRS. 

The last approach, a null-space method, is based on the RQ decomposition 

A=%R]Q, (12) 

where i? is a m-dimensional upper triangular matrix and Q is an n-dimensional orthogonal 
matrix. This RQ decomposition is obtained by using the LAPACK routine DQERQF. Premul- 
tiplying equation (6) by Q, we get QBQ^Qx + QA^y = Qb or 



^11 


B12 




Xl 


+ 
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B22 




X2 
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^2 



(13) 



where B = QBQ^, b = Qb and x = Qx. Matrices B and b are obtained by using the LAPACK 
routine DORMRQ. Now Rx2 = c by (7), (12) and Bnxi = bi-Bi2X2 by (13). The last equation 
is solved by using the LAPACK routine DSYSV. Finally R^y = 1)2 — B21X1 — B22X2 and x = Q^x 
can be obtained by the LAPACK routine DORMRQ. 

Notice that ABS algorithms were implemented in their basic form without partitioning into 
block or other special adjustements serving for speed increase as done in the LAPACK software. 

For each selected problem, the type of matrix and the dimension is given. Furthermore, both 
the solution and the residual errors together with the detected rank and the computational time 
are given for each tested algorithm. Computational experiments were performed on a Digital 
Unix Workstation in the double precision arithmetic (machine epsilon equal to about 10~^^). 

In the Appendix we give detailed results. The following tables give synthetic results for 
the 24 tested problems, the number at the intersection of the i-th row with the k-ih column 
indicating how many times the algorithm at the heading of the i-th. row gave a lower error than 
the algorithm at the heading of the A;-th row (in case there is a second number, this counts the 
number of cases when difference was less that one percent). 
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solution error - 



24 KKT linear systems 



mod. impl. impl. lu range null 

huang lu8 lu9 lapack space space total 
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impl . Iu9 
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lu lapack 


13 


14 


15 
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range space 
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null space 
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12 


13 
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20 
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residual error - 24 KKT linear systems 

mod . impl . impl . lu range null 

huang lu8 lu9 lapack space space total 

mod. huang 15/1 16 11 20/1 6 68/2 

impl.luS 8/1 15/2 8 19 5 55/3 

impl.lu9 8 7/2 7/2 20 6 48/4 

lu lapack 13 16 15/2 22 8/2 74/4 

range space 3/1 5 4 2 5 19/1 

null space 18 19 18 14/2 19 88/2 



Prom the above tables and the results in the Appendix we can state the following conclusions: 

(1) Implicit LU methods described in Subsection 2.3 are extremely suitable in term of total 
time cost for solving KKT systems with m ~ n. If m « n, then a direct solution of a 
complete KKT system by using the Bunch-Parlett decomposition is faster. 

(2) The method based on equation (9) is faster if m ~ n while the method based on equation 
(11) is more efficient if m « n. 

(3) The modified Huang method is usually very slow. But in the extremely ill-conditioned cases 
it gives the best accuracy. Moreover, the modified Huang method is able to detect the 
numerical rank correctly, which can lead to the substantial decrease of the computational 
time as can be observed from experiments with the matrix IDF2. 

(4) In term of accuracy the range space method is the least accurate; other methods have 
a comparable accuracy, with a marginal advantage for lu lapack (except for the very ill 
conditioned problems, where mod. huang gives the best results, up to four orders better). 
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4 Appendix: Test results for KKT linear systems 



matrix dimension method solution residual rank time 

n m error error 
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impl . Iu9 
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0. 


.18D- 


13 


1900 


19. 


,00 


IR500 


1000 


900 


lu lapack 
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range space 
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0, 


.52D- 


11 


1900 


79. 
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IR500 
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900 


null space 
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-06 


0, 


.20D- 
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1900 


85. 
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condition number: 
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mod . huang 


0. 


,32D- 


■07 


0. 


.90D- 
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1800 


188. 


,00 


IR500 
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600 


impl . Iu8 


0. 


,83D- 
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0, 
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impl . Iu9 
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lu lapack 


0. 


,49D- 


■07 


0. 
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condition number: 
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lu lapack 
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mod . huang 
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. 74D-I-00 


0, 


.75D- 


•11 


1900 


79. 


,00 


IR500R 


1000 


900 
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Test results for KKT linear systems - continued 
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mod . huang 
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lu lapack 
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Test results for KKT linear systems - continued 



matrix dimension method solution residual 

n m error error 
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Test results for KKT linear systems - continued 



matrix dimension method solution residual rank time 

n m error error 
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